
%% initialize
clear all;
close all;
rng(328120)

%% read in data

% choices

rawdata = csvread('../../01_data/05_temp/hungarian_input.csv',1,0);
VA1 = rawdata(:,1);
VA2 = rawdata(:,2);
N1  = rawdata(:,3);
N2  = rawdata(:,4);
exp = rawdata(:,5);
s   = rawdata(:,6);
N1c = rawdata(:,7);
N2c = rawdata(:,8);

I = length(VA1);

clearvars rawdata;

% allocations with switching costs


lambdastart = 0;
lambdaend = 25;
lambdainc = .05;
nlambda = (lambdaend-lambdastart)/lambdainc+1;

output = zeros(I,nlambda);

num_reassigned = zeros(nlambda,1);

for lambda=lambdastart:lambdainc:lambdaend
    
    lambda

n = round((lambda-lambdastart)/lambdainc+1);

C1 = zeros(I);


for i=1:I
    for j=1:I
        C1(i,j) = -(VA1(i)*N1(j) + VA2(i)*N2(j))-lambda*(i==j);
    end
end

[iind_switch,jind_switch] = linear_sum_assignment(C1);

tempreassigned = 0;
for j=1:I
   output(j,n) =  VA1(iind_switch(j))*N1(j) + VA2(iind_switch(j))*N2(j);
   tempreassigned = tempreassigned + (find(iind_switch==j)~=j);
end

num_reassigned(n) = sum(tempreassigned);


end

output_per = sum(output,1)./(sum(N1+N2)*ones(1,size(output,2)));

first_nonew = find(num_reassigned==0,1,'first');

frac_reassigned = num_reassigned/I;

plot(frac_reassigned(1:first_nonew),output_per(1:first_nonew))

switch_cost_mat = [output_per(1:first_nonew)',frac_reassigned(1:first_nonew)];

writematrix(switch_cost_mat,'../../01_data/05_temp/gains_switch_costs.csv','WriteMode','overwrite');

% best allocation

C1 = zeros(I);

for i=1:I
    for j=1:I
        C1(i,j) = -(VA1(i)*N1(j) + VA2(i)*N2(j));
    end
end

[iind_best,jind_best] = linear_sum_assignment(C1);

writematrix(iind_best,'../../01_data/05_temp/assignments_best.csv','WriteMode','overwrite');

% worst allocation
C1 = -C1;

[iind_worst,jind_worst] = linear_sum_assignment(C1);

writematrix(iind_worst,'../../01_data/05_temp/assignments_worst.csv','WriteMode','overwrite');

% best within experience allocation

C1 = zeros(I);

for i=1:I
    for j=1:I
        C1(i,j) = -(VA1(i)*N1(j) + VA2(i)*N2(j));
        if exp(i)~=exp(j)
           C1(i,j) = 100000; 
        end
    end
end

[iind_wiexp,jind_wiexp] = linear_sum_assignment(C1);

writematrix(iind_wiexp,'../../01_data/05_temp/assignments_wiexp.csv','WriteMode','overwrite');

% best within school allocation

C1 = zeros(I);

for i=1:I
    for j=1:I
        C1(i,j) = -(VA1(i)*N1c(j) + VA2(i)*N2c(j));
        if s(i)~=s(j)
           C1(i,j) = 100000; 
        end
    end
end

[iind_wis,jind_wis] = linear_sum_assignment(C1);

writematrix(iind_wis,'../../01_data/05_temp/assignments_wis.csv','WriteMode','overwrite');

% best allocation for m1

C1 = zeros(I);

for i=1:I
    for j=1:I
        C1(i,j) = -(VA1(i)*N1(j));
    end
end

[iind_maxm1,jind_maxm1] = linear_sum_assignment(C1);

writematrix(iind_maxm1,'../../01_data/05_temp/assignments_maxm1.csv','WriteMode','overwrite');

% best allocation for m2

C1 = zeros(I);

for i=1:I
    for j=1:I
        C1(i,j) = -(VA2(i)*N2(j));
    end
end

[iind_maxm2,jind_maxm2] = linear_sum_assignment(C1);

writematrix(iind_maxm2,'../../01_data/05_temp/assignments_maxm2.csv','WriteMode','overwrite');

% actual allocation

iind_actual = (1:I)';

writematrix(iind_actual,'../../01_data/05_temp/assignments_actual.csv','WriteMode','overwrite');

% random allocation

S = 100;
iind_random = zeros(I,S);

for s=1:S
    random_sort = rand(I,1);
    [~,iind_random(:,s)] = sort(random_sort);
end

writematrix(iind_random,'../../01_data/05_temp/assignments_random.csv','WriteMode','overwrite');


% best allocation, with classroom variation

C1 = zeros(I);

for i=1:I
    for j=1:I
        C1(i,j) = -(VA1(i)*N1c(j) + VA2(i)*N2c(j));
    end
end

[iind_best_NC,jind_best_NC] = linear_sum_assignment(C1);

writematrix(iind_best_NC,'../../01_data/05_temp/assignments_best_NC.csv','WriteMode','overwrite');

% save long data


output = zeros(I);
itidx = zeros(I);
jtidx = zeros(I);
for i=1:I
    for j=1:I
        output(i,j) = (VA1(i)*N1(j) + VA2(i)*N2(j));
        itidx(i,j) = i;
        jtidx(i,j) = j;
    end
end

output = output(:);
itidx = itidx(:);
jtidx = jtidx(:);

save('../../01_data/05_temp/full_sample_output', 'output','itidx','jtidx','N1','N2','VA1','VA2');
